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Abstract 

We calculate the two-photon ionization of clusters for photon energies near 
the surface plasmon resonance. The results are expressed in terms of the ion- 
ization rate of a double plasmon excitation, which is calculated perturbatively. 
For the conditions of the experiment by Schlipper et al., H, we find an ion- 
ization rate of the order of 0.05-0.10 fs _1 . This rate is used to determine the 
ionization probability in an external field in terms of the number of photons 
absorbed and the duration of the field. The probability also depends on the 
damping rate of the surface plasmon. Agreement with experiment can only 
be achieved if the plasmon damping is considerably smaller than its observed 
width in the room-temperature single-photon absorption spectrum. 

I. INTRODUCTION 

The electromagnetic response of alkali metal clusters shows a very strong surface plasmon 
resonance PJ] , but the interactions of the plasmon with other degrees of freedom are not well 
understood. One interesting question is the nonlinearities associated with multiple plasmon 
excitations-how weakly do they interact with each other? Some physical processes can 
be sensitive to nonlinearities; for example ionization may be energetically impossible for 
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individual plasmons but allowed for states with multiple plasmon excitations and therefore, 
ionization rates may depend on the degree of nonlinearity. Recently an experiment was 
reported observing the ionization probability with different photon field durations ||. The 
photon energy was such that ionization is energetically possible only if at least two photons 
are absorbed. In this work we ask whether the observed ionization can be interpreted as 
electron emission by a two-plasmon state within a simple theory based on the jellium model 
of the electronic structure. We will use standard many-body perturbation theory, describing 
the surface plasmon in RPA as a particle-hole excitation. The plasmon description and 
details of the jellium model are given in Sect. 2 below. 

Our calculation may be viewed as a semianalytic approximation to time-dependent local 
density approximation (TDLDA), in which the electron dynamics is treated entirely by a 
common potential field. The TDLDA has been well developed for the high-field response of 
atoms 0-§f, and is now being applied to sodium clusters ||. Unfortunately, the full TDLDA 
is computationally difficult and rather opaque, in contrast to the perturbative approach that 
allows important quantities to be calculated directly. 

From the point of view of surface plasmon dynamics, a very important quantity is the 
ionization rate of a two-plasmon excited state. Haberland et al. || interpreted their mea- 
surements under the assumption that this rate is fast on a time scale of 10 fs, and we wish 
to see whether that can be justified theoretically. The two-plasmon ionization rate is cal- 
culated in Sect. 3. However, the ionization can take place without the mediation of the 
plasmon. Also, the plasmon can be excited as a virtual state in which case the connection to 
the two-plasmon decay formula is unclear. We present in Sect. 4 a more general treatment 
of the ionization process that includes these effects and allows the role of the plasmon to 
be isolated from other aspects of the ionization process. The important role of plasmon in 
screening and enhancement of the external field is made explicit in the formulas discussed 
there. 
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II. THE ELECTRONIC STRUCTURE MODEL 



In this section we discuss the Hamiltonian model and the treatment of the surface plas- 
mon. We will need single-electron wave functions and energies, which we calculate as follows. 
We first obtained the solution of the self-consistent jellium model using the computer code 
"jellyrpa" M. The jellium background charge is assumed to be uniform in a sphere of radius 
R = TsN 1 ^. Here r s = 3.93 a.u. corresponds to density of charge equal to the bulk density 
of atoms in sodium metal. It is more convenient for recreating the wave functions to work 
with analytic models of the potential, so we fit the self-consistent jellium potential to a 
Woods-Saxon shape. Specifically, we take the electron potential to be 

where V c {r) is a Coulomb field associated with the positive charge distributed uniformly in 
the jellium sphere, 

e 2 

VJr) = — , r > N 1/3 r s 
r 

The parameters that fit this potential to the self-consistent one are Vq = 5.71 eV, Rq = 
10.548 A, and a = 0.635 A. The occupied energy levels of this potential are within 0.2 eV 
of the self-consistent potential, which is certainly within the accuracy of the jellium model. 
We find that the cluster has an ionization potential of 4.5 eV. Under the conditions of the 
experiment [2] using photons of 3.1 eV, two photons are required for ionization on energetic 
ground. The single-electron spectrum is shown in Fig. 1. We use these orbitals and energies 
in the RPA and ionization calculations. 

The RPA surface plasmon might also be calculated numerically with the code "jellyrpa" , 
but in the interests of developing analytic formulas we adopted a more schematic approach. 
We take the interaction between electrons to have a separable form f|, 



v r, r 



= (r) • f (r) , (3) 



where f is a three-dimensional vector with components f) t ( r ) = /M^TCO- Then the energies 
of the RPA excitations satisfy the dispersion relation 

where e is a single-particle energy and p, h label particle and hole orbitals. Due to the 
spherical symmetry, the solutions uj n of the dispersion relation are independent of fi. The 
matrix element (n/ilfjO) between the ground state and a one-phonon state of energy u n is 
given by 

We shall particularly require the transition potential v nfJL associated with the creation of the 
plasmon. This is given by 

v ni i = Kf^(r)(n/i|f M |0) . (6) 

In the spherical jellium model, the surface plasmon can be roughly described taking the 
interaction of dipole-dipole form. For an excitation along the z-axis, the field is 

f (r) =z. (7) 

Assuming that the transition density of the plasmon is concentrated at the surface at radius 
R, the strength of the interaction is obtained from the multipole expansion of the Coulomb 
interaction as 

The dispersion relation can then be solved analytically H] in the limit uo » (e p — e^) making 
use of the TRK sum rule. The result is the simple Mie surface plasmon formula, 

The resulting energy is about 25% higher than the empirical value for sodium clusters, 
uj m 2.75 eV. The RPA can be made to fit this value for iV = 93 by renormalizing the 



coupling strength by k = 0.52k c . However, the transition potential of eq.(|J) calculated with 
eq.(|7|) is linear in r whereas TDLDA calculations without separable assumption do not yield 
this behavior, as shown in Fig. 2. A simple improvement over the linear form eq. ([?[) is 
the dipole field associated with a charge distribution localized on the surface of the jellium 
sphere [110] . A surface charge produces a radial field of the form 



/(r) = r r < R 



= ^ r>R. (10) 

This is plotted as the dashed line in Fig. 2; it clearly is much closer to the actual TDLDA 
transition potential. With this choice, the empirical position of the resonance is obtained 
by using the coupling k = 0.6k c . This is very close to the previous one, showing that the 
modified form factor has only a small influence on the plasmon properties. We will see that 
it is much more important in the ionization process. 

In the experimental photoabsorption spectrum, the plasmon has a width of about 0.5 eV. 
This finite width requires theory beyond RPA, which produces only zero-width excitations 
below the ionization threshold. Since it is not easy to incorporate the other degrees of 
freedom that are responsible for the width, we will treat the width phenomenologically. In 
discussing the response in general, it is useful to consider the dynamic polarizability a(u). 
This is given by 

x . 2e 2 (n\z\0) 2 u n 

a { u} ) — / — ^ 7 ^7^7 (11) 

where n labels the true excitations and 5 is a small quantity. A simple prescription is to 
take a single pole for the plasmon, taking into account the finite width by the replacement 
5 — > T n /2 for a width T n . 

The imaginary part of a is related directly to the photoabsorption cross section a by the 
formula 

a = 4vr^Im a(u). (12) 



Given Im a (a;), the real part can then be computed from the Kramers-Kronig relation, 

2u' 



If 00 , ( 2u \ 

Rea(u;) = — du'lmaiu') — . (13) 

n Jo \uj 2 — uj' 1 I 



Applying this to the experimental data of ref . [O , we find the imaginary part of a graphed as 



the solid line in Fig. 3. This is compared with the single-pole approximation with parameters 
uj n = 2.75 eV and T n = 0.5 eV. A modification of the jellium model was proposed in ref. 



12j introducing a soft-edged surface in the distribution of the background charge. We 
also calculated the full RPA response for soft jellium model, calculated using the program 
"jellyrpa". In this case the empirical width can be reproduced with a smaller external 
width parameter. The dashed line shows the fit with T n = 0.3 eV. Both these models give 
a reasonable but not quantitative description of the data. The soft jellium model has the 
advantage that the plasmon can be moved to lower frequency without adjusting the coupling 
strength. However, it predicts too low an ionization potential, which makes it unsuitable 
for the autoionization calculation. 

The corresponding comparison for the real part of a is shown in Fig. (4). Here the theory 
is quite robust, and we can rather confidently estimate Re a at the energy of interest (3.1 
eV) to be about 3000 A 3 . 

III. TWO-PLASMON AUTOIONIZATION RATE 

The many-body perturbative graphs for M p h, the interaction matrix element between 
the two-plasmon excitation of the ra/i mode and the final configuration with a hole h and 
the electron in a continuum state p , is shown in Fig. 5. The labels p', h! stand for particle 
and hole states, respectively. Algebraically, the matrix element is given by 

M ph = s/2 



I ( h \ V nn\p')(p'\Vnn\p) \ _ f( h \ V nfl\h'){h'\v nfl \p) 



— V + C?i / h t \ — £p + £h 



w n — U' + £h 

where t>„ M is defined in eq. (^). The factor \/2 accounts for the statistics of the two-plasmon 
initial state. The two graphs can be combined in one sum over both particles and holes as 
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shown in the second line, making use of the fact that the matrix element is only required 
on shell, i.e. with e p — = 2u n . The primed indices %' indicate particle or hole orbitals, 
depending on the direction of the arrow. The ionization width T e = hw e , where w e is the 
ionization rate, is given by the Golden Rule formula, 



,dn r 



T e = 2ttJ2 \M ph \ 2 ^£5{2uj n - e p + e h ) 



ph 



4™ 4 ( n|f M |0> 4 £ \K ph \ 2 ^6(2uj n - e p 4 r,,) 



ph 



dE 



(15) 



Here dn p / dE is the density of states of the continuum electron. We have also separated out 
the excitation amplitude for a field f , 



ph 



E 



~ £%> + Ch 



(16) 



The sums in eq. (|15|) can be reduced in size by making use of the angular momentum 
symmetry of the orbitals. Labeling the angular momentum quantum numbers / and m, we 
may express the m-dependence of the matrix elements as 

1 I, 



(p',m p/ \^\p, m p ) 



-l)V-"v 



m p i /i m p 



(p'\\f\\p) 



(17) 



where the reduced matrix element (o||f||6) is defined as 

la 1 lb \ /-00 



(a||f||6) 



(2l a + l)(2l b + 1) 



f(r)ip a (r)(p b (r)r 2 dr 



(18) 



/ Jo 

in terms of the radial wave functions </?j. The sum over magnetic quantum numbers m p ^ 
implicit in eq. (|l^) can be evaluated in terms of a 9-j symbol |l4j] in which the total angular 
momentum L carried by the two photons appears. The result is 

E I^| 2 = 2E-J37 1 1 

all m 



Ej -e h )u- (tj - e h ) 

f 1 1 L \ 



(fe||f||i>01|f||p><*||f||fe>(p||f||*> 



1 1 V 

x E L 1 " 

L=0,2 VO 



3 p ^ 



(19) 



The factor of 2 arises from the two-fold spin degeneracy of the occupied orbitals. In carrying 
out the calculation one also has to fix the normalization of the continuum radial wave 
function. A convenient choice is rip p — > sin(kr + 8) at large r, giving dn p /dE = 2m/(iikh 2 ). 

The results for the autoionization of Na^ are given in Table I. We chose the plasmon 
parameter k in two different ways, requiring the plasmon resonance energy uj n to be either at 
the experimental position of 2.75 eV, or at the energy of the absorbed photons, 3.1 eV. The 
particle-hole states taken into account in the calculations include electron jumps up to three 
harmonic oscillator shells. We first discuss the results in the case of undamped excitations 
(5 = in Table I). The upper half of Table I shows the values obtained using the linear 
dipole field f = z (eq. (J?])). The resulting widths T e are extremely small for both choices 
of u n , and they would be hard to reconcile with experiment. This led us to reexamine our 
simplifying assumption about the shape of the separable particle-hole interaction. Since the 
choice of eq. fllCf ) gives a better transition potential (see fig. 2) we use it from now on instead 
of eq. (|7|). As shown in the lower half of Table I, the resulting widths are much larger and 
they seem to give a possibility to explain the data. Indeed, they correspond to ionization 
times of the order of 5.5 fs to 7.5 fs which is comparable to the estimated plasmon lifetime 
of 10 fs §. 

However, the calculated results cannot be considered reliable because they are quite 
sensitive to the single-particle energies involved in the transition. In eq. ([16]) several states 
i' give quite small energy denominators (see Fig. 1) and therefore they yield abnormally 
large contributions. However, it is not consistent to neglect the damping of the excitations 
in the perturbative calculation when the energy denominators are small. As we did with 
the plasmon in Sect. 2, we here add a finite imaginary term 5 to the energy denominators of 
eq. (|16|). In Fig. 6 we show the dependence of \Kph\ 2 on 5 which is seen to be moderate 
for 5 in the range 0.1 - 0.2 eV. We see that the damping reduces the ionization width by 
roughly a factor of two when 5 = 0.1 eV which seems a reasonable value. Although the 
dependence on the transition field and on the damping makes the calculation uncertain by 
a factor of two or so, we at least see that the result is rather insensitive to the particular 
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model for the plasmon excitation energy. This is nice in that it means that the theoretical 
uncertainty here is not a hindrance to doing the calculation. 

The results correspond to a ionization lifetime r e of the two-plasmon state in the range of 
5 to 10 fs. This is of the same order of magnitude as the plasmon lifetime which is estimated 
to be about 10 fs |2|]. In the present model the ionization process is not very fast contrarily 
to the assumption made in the interpretation of ref. ||, but it is fast enough to allow for 
the ionization process to compete with the plasmon damping. 



IV. IONIZATION IN A LASER FIELD 

In this section we wish to apply the previous results to ionization in a laser field. Thus, 
we consider the ionization as a multistep process, in which the photons are first absorbed 
to make plasmons and then the plasmons interact to eject an electron. A simple physical 
argument can be made to obtain a formula for the ionization, which we will then justify 
more formally. Let us define the absorption rate for photons R v and the damping width 
for the plasmon T n . In the steady state the balance between the creation and absorption of 
plasmons gives a mean number of plasmons n satisfying 

n = ^ (20) 

Taking the distribution of numbers as Poissonian, the mean number of pairs is then n 2 /2. 
The ionization rate R e is related to the two-plasmon ionization width T e by 

n 2 

Re = — W P , 



2 

= |#r e , (2i) 

where T e and w e have been introduced in the preceding section. 

A more formal derivation of this formula may be made from the graph of Fig. 5 as 
follows. We add to the graph matrix elements of the external field between the zero- and 
one-plasmon states, 



V n = (n\V ext \0) . (22) 

The plasmon propagator itself can be approximated by + u n + iT n /2) for u close to 

uj n . Then the ionization rate for the graph is given by 

2 



pa 



-co + oo n + iT n /2 



' dn p 
~dE 



i v r 

ft (K - + (r n /2) 2 ) 2 2 ' 1 J 

On the other hand, the photon absorption rate can be calculated as the imaginary part of 
the self-energy associated with the coupling V n , 

Ru = ^K 2 Im- ' 



h n —uj + uj n + iT n /2 
1 V 2 T 

= 1^ . (24) 

Eq. ( pi]) may now be obtained by combining the last two equations. 

For comparing with experiment, it is convenient to express eq. ( PH) in terms of the number 
of ionizations per cluster N e = R e T and the number of photons absorbed N v = R U T, where 
T is the time duration of the laser pulse. We obtain for iV e 

= ^ • (25) 

n 

The experiment [2] observed not only ionization but considerable evaporation of atoms 
from the clusters. We note that atomic motion takes place on a much longer time scale 
than electronic motion. Most of the evaporated particles are emitted from the cluster after 
it has reached thermal equilibrium; the statistical evaporation theory gives lifetimes in the 
nanosecond regime for the conditions of the experiment [TJJ. This is longer than the fast 
laser pulse duration (T = 140 fs), so we may think that the Nag" 3 cluster remains whole 
before being ionized. The first line of eq. (25) may be estimated from Fig. 2 of ref. M. The 
broad peak in the middle of the mass spectrum represents ionized clusters and has roughly 
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1/3 of the area of the Nag~ 3 peak on the right; thus N e m 1/3. The number of photons 
absorbed is given as N v » 6. Thus we estimate^ P e « 0.014. 

To evaluate the second line of eq. (25) we have the calculated T e from Sect. 3 and the 
experimental pulse duration T, but we still do not have T n , the damping width of the 
plasmon. As discussed in Sect. 2, this quantity is beyond the scope of RPA, but an upper 
bound to r„ is given by the empirical photoabsorption spectrum. Taking the values ui n =2.75 
eV and T n = 0.5 eV the single-pole fit to the plasmon and 5 = 0.2 eV in the perturbative 
energy denominators, we find P e = 1.0 x 10~ 3 , which is much too small. On the other 
hand, the experimentally quoted lifetime of the plasmon is 10 fs which corresponds to a 
considerably smaller value of T n . Adopting this value of T n gives P e = 5.8 x 10 -2 . The 
width T e depends on the particular cluster considered since the energetics of the emitted 
electrons will change somewhat for different clusters. 



V. GENERAL THEORY 

The theory in the last section assumed that the energy transfer to the electron was 
indirect, first producing plasmons which then autoionize. In fact the photon could be ab- 
sorbed directly on the electron. The direct absorption is implicit in the TDLDA, and can be 



taken into account as well in the perturbative theory of photon absorption ||16|| . The general 
expression for a second-order transition from a state i to / is given by 

^fe W '§. (26) 

The second order matrix element Kfi is similar to eq. (16) with the linear field (7), except 

the particle orbitals are replaced by many-body states f: 

= y, (f\z\i')(i'\z\i) ^ 
., lo — Ei> + Ei 



1 Fitting more extensive data with a rate equation, the authors of ref. Q find a somewhat larger 
probability, P e ~ 0.06 
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In eq. (EH) £ is the amplitude of the time- dependent electric dipole field, £(t) = EtfL{e~ lu)t + 
e luJt ). It is related to the laser intensity by Sq = 2ttIuj/c. Taking the many-body states 
as the simple particle- hole configurations, eq. ( p7|) reduces to eq. (|l6l). By itself this would 
only give the ionization probability of the electron in the external field, i.e. without plasmon 
effects. The many-body physics is included by considering higher-order perturbations in the 
wave functions including excitation of the other electrons by the active electron. The result 
is to replace the external field So in eq. fl26|) by an effective field such that 

eS eff z = eS z -ej v(r,r')U(r' \r n )£ z n d 3 r'd 3 r n , (28) 

where II is the response function of the cluster. 

To make the connection with the previous approach, we approximate the response by a 
single pole and use the separable approximation. The effective field becomes 

eS eff z = eS z - e£ O Kz(n\z\0)(0\z\n) ( 1 + — 1 , ) 

= e£ z(l - ^a(uj) . (29) 



Eq. (^) can now be obtained from eq. ( p9|) by dropping the external field contribution as well 
as the nonresonant term in the polarizability. In the second line, we express the polarization 
effects as a multiplicative factor, which can be interpreted as an effective charge 

e eff = e(l - -^a(u)) . (30) 

We may use this expression to assess the relative importance of the external and induced 
fields. First note that there is a complete cancelation of the two terms in eq. (30), implying 
complete screening, if we take k = k c and a = R 3 . This last is just the classical polarizability 
of a conducting sphere, and it may also be derived putting the Mie resonance eq. (0) into 
the polarizability formula eq. (JTTJ) . 

For the present purposes we assume in eq. (29) u n = 2.75 eV and k = 0.52k c as in Sec. 2. 
We also take for T n a width corresponding to a lifetime of 10 fs. The result at u = 3.1 eV is 

4«(3.1 eV) = -3.3 + i0.3 (31) 



e 2 
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Thus the plasmon enhances the field approximately by a factor of 4.3, showing that the 
induced field is indeed dominant. This result is insensitive to T n if it is small, but would 
decrease if the width were as large as the measured width of the optical absorption peak. In 
a study of Nag" in the TDLDA, ref. [T7[ obtained enhancement factors on resonance in the 
range 5-8, which is the same order of magnitude. 

The effect on the ionization rate goes as the fourth power of the field, 

K,a(3.1 eV) ./ , , 

1 ^—5 -| 4 = 320 , (32) 

e 2 

Finally , we use this result to make an improved calculation of the ionization probability P e 
introduced in eq. (|2"5]) . The number of emitted electrons is calculated as N e = Wf(T where 
Wfi is given by eq. ( p6|) with e replaced by e e //. The number of of pairs of absorbed photons 
is given by N pair = (IaT) 2 /2, a being the photo- absorption cross section. Computing a from 
eq. flT2D we find P e = 1.25 x 10 -2 , a magnitude comparable to that found in the preceding 
section. However, the cross section corresponding to the polarizability in eq. (|3~TD is much 
too small at 3.1 eV, as is clear from the width discussion in Sect. 2. On the other hand, 
if we take the width parameter from empirical single-pole fit, T n = 0.5 eV, the probability 
comes out very small, as discussed earlier. 



VI. CONCLUSION 

In this work, we have derived the theory of cluster ionization by multiple photons of fre- 
quency near that of the surface plasmon. The weak coupling between the surface plasmons is 
the driving interaction for the two-photon ionization process, and a perturbative framework 
with respect to the ionized electron seems reasonable. The plasmon-induced mechanism 
can be derived from the general perturbative formula using the higher order contributions 
associated with the screened interaction. Unfortunately, the formula depends quadratically 
on the damping rate of the plasmon, which is still not fully understood. The rates obtained 
for Na^ are of the order of tens of femtoseconds, which is the same time scale as other 
relaxation processes. 
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We have used the jellium model in the calculations, and it is unclear how realistic the 
model is. We found that the interaction must be treated more accurately than in the small- 
r separable approximation, but we have not examined the most sophisticated treatment of 
the interaction which would include the exchange interaction without approximation. The 



Landau damping of the Mie resonance is much larger with more realistic interactions [ITS 
A major problem of the jellium model is that the damping is too small. Ionic scattering 
(called "interband transitions" in condensed matter physics) would increase the spreading 
width, and the lack of ionic scattering in the jellium model is in general a serious deficiency. 
However, in the case of Na clusters the realistic Hamiltonian gives a very similar spectrum to 



the jellium model |19[], giving some credibility to the model. Nevertheless, it would be inter- 
esting to see what the effects of the ionic scattering are in the second-order ionization. For 
small Na clusters, more realistic calculations are becoming available of high field ionization 



using the TDLDA [20], and it would be interesting to compare. 
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TABLES 

TABLE I. Two-plasmon ionization widths in Na^. The upper half and lower half of the table 
correspond to results calculated with the coupling (0) and ([jl]), respectively. The quantities r e are 
defined in the text. 
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FIGURES 
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FIG. 1. Single-particle levels in the jellium model of Na^. The Fermi level is indicated with 
a dotted line. The arrows show a two-step transition with a particularly small matrix element in 
eq. (16). 
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FIG. 3. Imaginary part of the dynamic polarizability of Na^: empirical from ref. and 
eq. ( li~2| ) (solid line); single-pole approximation (dashed line); RPA with soft jellium model (dotted 
line). 
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FIG. 4. Real part of the dynamic polarizability for Na^Jg. Empirical from Kramers-Kronig 
relation eq. (|i~3|), solid line; single-pole approximation (dashed line); RPA with soft jellium model 
(dotted line). 
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FIG. 5. Perturbation theory graphs for second-order ionization. 
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FIG. 6. Second-order term ne^ /2\K^\ 2 dn f / dE as a function of 5. 
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